Capítulo 2: Sobre a natureza das relações

Relações lineares

Uma relação linear entre duas variáveis indica que elas estão correlacionadas em uma proporção constante para qualquer intervalo.

Isto é, valores maiores de massa correspondem a um aumento proporcional em energia. O valor de \(c^{2}\) expressa essa proporção constante.

Exemplo: uma molécula de água pesa aproximadamente \(m_{H_{2}O} =2.992×10^{-23} g\). Portanto, a energia associada é \(E_{H_{2}O} = 2.992×10^{-23}*8.988*10^{16} \sim 2.689^{-6}J\). Se triplicarmos o número de moléculas de água, o mesmo acontecerá com a energia associada: \(E_{3H_{2}O} = 3*E_{H_{2}O}\).

Se a correlação é positiva, incrementos em \(x\) serão proporcionais a incrementos em \(y\). Se a correlação é negativa, incrementos em \(x\) serão proporcionais a decréscimos em \(y\).

Num cenário perfeito, se sabemos que há uma relação linear entre variáveis, precisamos de apenas duas observações para descobrir proporção entre elas. Esse problema é idêntico ao de encontrar a inclinaçaõ da reta que passa por dois pontos. É de fácil resolução usando técnicas elementares.

    >library(ggplot2)
    >ggplot()+
    geom_point(mapping=aes(x=1,y=2))+
    geom_point(mapping=aes(x=2,y=4))+
    xlim(0,3)+ylim(0,5)

\(y=\beta*x\)
\(a=(1,2); b=(2,4) \rightarrow \beta = 2\)

    >ggplot()+
    geom_point(mapping=aes(x=1,y=2))+
    geom_point(mapping=aes(x=2,y=4))+
    xlim(0,3)+ylim(0,5)+
        geom_abline(slope = 2)

Erros e aleatoriedade

Controlando fatores experimentais, as relações descritas são bastante precisas. Em um cenário sem atrito com superfícies e com o ar, os erros de medida obtidos com \(F=ma\) são muito baixos.
Entretanto, nem sempre isso é verdadeiro.
Primeiro, podemos sofrer interferência de variáveis desconhecidas.

Imaginemos um conjunto de medidas antropométricas, com altura e peso e indivíduos.
É esperado que a altura de um ser humano esteja relacionada com seu peso. Entretanto, outras características não medidas, como percentual de gordura total, podem interferir nos valores finais. Normalmente, tratamos essas flutuações como erros aleatórios¹.

Podemos simular este cenário partindo de variáveis idênticas e adicionando ruído aleatório.

    >set.seed(2600)
    >a <- seq(1:100)+rnorm(n=100,sd=3)
    >b <- seq(1:100)+rnorm(n=100,sd=3)

    >cor_data <- data.frame(a,b)
    
    >ggplot(cor_data,aes(x=a,y=b))+
    geom_point()

O resultado sugere que há uma forte relação linear entre \(x\) e \(y\). Por outro lado, notamos que é impossível para uma reta cruzar todos os pontos. A seguir, vamos investigar como quantificar a correlação linear, assim como encontrar a reta que minimiza a distancia para todas as observações.

Com essas ferramentas, podemos estender nossas inferências. Além de comparações, teremos noções sobre a magnitude de uma relação, assim como poderemos prever o valor esperado para novas observações.

O coeficiente de correlação produto-momento de Pearson, ou, simplesmente, \(\rho\) de Pearson.

O coeficiente de correlação \((\rho)\) de Pearson é um número real garantidamente² entre -1 e 1. Expressa a magnitude e o sentido de uma relação linear, sendo -1 uma relação inversa perfeita e 1 uma relação direta perfeita.

Para os dados que geramos, a correlação é quase perfeita: \(\rho = 0.989\).
O coeficiente possui produto-momento em seu nome, pois usa uma abstração originalmente empregada na física: o momento.

Dê-me um ponto de apoio e eu moverei a Terra

Um breve mergulho na física: Momentos

³

Para adquirir uma intuição sobre o coeficiente, é interessante resgatar o conceito físico, originalmente descoberto por Arquimedes. Embora não tenha inventado a alavanca, ele descreveu os princípios matemáticos por trás dela.

Em Sobre o equilíbrio dos planos, Arquimedes declara que Magnitudes ficam em equilíbrio quando em distância reciprocamente proporcional aos seus pesos.

Essa é a conhecida Lei da alavanca. Dado um ponto de apoio e um plano sobre ele, aplicamos uma força em qualquer local do plano. O momento (torque) resultante é o resultado da multiplicação da grandeza física \((F)\) pela distância até o ponto fixo \((d)\).
\(M = F*d\)

Supondo uma força constante, quanto mais nos afastamos do ponto fixo, maior o momento resultante. Posteriormente, os físicos estenderam o conceito para outros domínios. Por exemplo, um objeto com cargas opostas \(-q\) e \(+q\) separados por uma distância \(d\) possui momento (momento dipolar elétrico) análogo: \(M = q*d\)
De uma maneira geral, falamos em momento ao multiplicarmos uma grandeza física por uma distância.

Momento resultante

Como o brinquedo acima fica equilibrado sobre apenas um ponto?

No caso da alavanca, vimos que cada força aplicada sobre o objeto está associada a um momento(torque). Sabemos que a gravidade atua sobre cada pedaço com massa compondo o todo. Podemos então calcular o momento resultante somando os momentos de todos os \(N\) pontos. Seja \(F_{i}\) a função descrevendo a força em cada i-ésimo: \[M = \sum_{i=1}^{N} F_{i}d_{i}\]

Um sistema, como o pássaro apoiado sobre o dedo, está em equilíbrio quando a soma dos momentos em relação ao ponto fixo é zero. Para cargas elétricas, o sistema é apolar quando o momento é zero. Na figura abaixo, vemos como a molécula de \(CO_{2}\) é apolar, enquanto a molécula de água é polar:

Os momentos descritos acima são expressões do primeiro momento, uma vez que a grandeza é multiplicada pela distância com expoente 1: \(d = d^{1}\).

Podemos calcular outros momentos, exponenciando o componente espacial (distância). Vamos estudar agora momentos de massa de um objeto unidimensional:
O momento zero de massa para um objeto é \(M_{0} = \sum_{i=1}^{N} m_{i}d_{i}^{0}\). Como \(d^{0}=1\), temos \(M_{0} = \sum_{i=1}^{N} m_{i}\), que é simplesmente a soma das massas de todos os pontos. O momento zero é a massa total.
O primeiro momento de massa para um objeto é \(M_{1} = \sum_{i=1}^{N} m_{i}d_{i}^{1}\) e determina o centro de massa em relação à dimensão \(d\).
O segundo momento de massa é \(M_{2} = \sum_{i=1}^{N} m_{i}d_{i}^{2}\) e é o momento de inércia. Corresponde à resistência do sistema a rotações⁴. O n-ésimo momento \(M_{n} = \sum_{i=1}^{N} m_{i}d_{i}^{2}\)

Generalizando momentos

Podemos generalizar ainda mais a abstração e calcular momentos de entidades abstratas, como variáveis aleatórias. Melhor: já fizemos isso anteriormente!

Seja \(f(x)\) a função que descreve uma distribuição de probabilidades para a variável,
Assim como o momento zero representa a massa total, aqui ele representa a probabilidade total (1).
O primeiro momento, correspondente ao centro de massa na mecânica estática, é a média.
O segundo momento, correspondente ao momento inercial, é a variância.
Os momentos terceiro e quarto normalizados informam sobre assimetrias (skewness) e peso de valores extremos (kurtosis).

Formalmente, seja \(d(x,x_{0})\) o valor da distância ao centro \(x_{0}\) de referência \((x-x_{0})\), o n-ésimo momento \(\mu_{n}\) é definido por:

\[\mu_{n}= \int_{-\infty}^{\infty} d(x,x_{0})^{n}f(x)dx\]

A integral acima corresponde à versão contínua da soma de partes discretas apresentadas antes para uma grandeza física, como a massa: \(M_{n} = \sum_{i=1}^{N}d_{i}^{n}m_{i}\)

Momento zero: \[\mu_{0}= \int_{-\infty}^{\infty} d(x,x_{0})^{0}f(x)dx = \int_{-\infty}^{\infty} f(x)dx = 1\].
Primeiro momento: \[\mu_{1}= \int_{-\infty}^{\infty} d(x,x_{0})f(x)dx\] , supondo centro em 0 \((x_{0}=0)\), temos a média, \[ \mu_{1} = \int_{-\infty}^{\infty} xf(x)dx\], também chamado valor esperado \(E[X]\).
Segundo momento: \[\mu_{2}= \int_{-\infty}^{\infty} d(x,x_{0})^{2}f(x)dx\]. Como vimos anteriormente, a soma dos quadrados dos desvios, nossa variância, \(\sigma^{2} = E[(x-\mu)^{2}]\).

Com os conceitos adquiridos em mãos, é fácil entender o \(\rho\) de Pearson.

Calculando correlações lineares

A noção de distância ou desvio se repetiu muitas vezes.
De fato, o coeficiente de correlação linear nasceu quando Francis Galton (1888) estudava numericamente dois problemas aparentemente distintos em antropometria ⁵:

1. Antropologia: Se recuperássemos de um túmulo antigo apenas um osso da coxa (fêmur) de um indivíduo, o que podereríamos dizer sobre sua altura?
2. Ciência forense: Com o intuito de identificar criminosos, o que pode ser dito sobre medidas diferentes de uma mesma pessoa?

Galton percebeu que, na verdade, estava lidando com o mesmo problema. Dadas medidas pareadas, \((x_{i},x_{i}')\), o que o desvio de \(x_{i}\) informa sobre o desvio de \(x_{i}'\)?

O fêmur recuperado do esqueleto de um faraó é 5 cm maior que a média. Quão distante da média esperamos que seja sua altura? Ingenuamente, podemos pensar que se uma das medidas é 1 desvio-padrão maior que a média, a outra também será 1 desvio-padrão maior. Galton percebeu que havia um armadilha nesse pensamento.

Apesar de haver uma relação entre as medidas, há também flutuações aleatórias: parte do desvio é resultante disso. Precisamos entender o grau de correlação pra fazer um bom palpite.

Então, propôs um coeficiente mensurando a relação entre desvios de variáveis. Se tamanho do fêmur e altura estão muito relacionadas, um fêmur grande sugere indivíduo igualmente alto. Caso contrário (baixa correlação), um fêmur grande (desvio alto) não implica grande estatura.

Para quantificar a relação, multiplicamos os desvios de cada par de medidas:
\[\sum_{i=1}^{N}(x_{i}-\mu_{x})(x_{i}'-\mu_{x'})\] A expressão acima expressa a covariância entre \(X\) e \(X'\) e será útil em outros contextos. A expressão lembra o cálculo do primeiro momento, porém cada desvio é multiplicado pelo desvio correspondente da medida pareada. Daí o nome coeficiente de correlação produto-momento.

Note que, se ambos os desvios concordam em sentido (sinal), o resultado da multiplicaçao será positivo. Pares consistentemente concordantes aumentam o valor da soma final. Se ambos os desvios discordam em sentido (sinal), o resultado será negativo. Pares consistentemente discordantes diminuem o valor da soma final.

Assim, podemos ter variáveis altamente correlacionadas positivamente negativamente, desde que o sentido da associação seja constante. Em contrapartida, se as medidas são ora discordantes e ora concordantes, os valores tendem a se anular na soma e o resultado se aproxima de zero.

Observar apenas a covariância é perigoso, pois os valores dependem da unidade de medida e da dispersão dos dados.

Calculamos o coeficiente de correlação de Pearson, normalizando⁶ a covariância ao dividí-la pelo produto dos desvios-padrão: \[\rho_{XX'}= \frac{cov(X,X')}{\sigma_{X}\sigma_{X'}}\] De forma extensa:
\[\rho_{XX'}= \frac{\sum_{i=1}^{N}(x_{i}-\mu_{x})(x_{i}'-\mu_{x'})}{\sqrt{\sum_{i}^{N}(x_{i}-\mu_{x})^{2}}\sqrt{\sum_{i}^{N}(x_{i}'-\mu_{x'})^{2}}}\] Uma boa notícia: \(\rho\) segue uma distribuição conhecida, a distribuição t, com n-2 graus de liberdade. Podemos usar as ferramentas anteriores para testar hipóteses.

Exemplo prático

O exemplo a seguir foi um feliz achado. Na época, o governo brasileiro discutia a necessidade da ampliar número de médicos para melhorar a assistência à saúde. Alguns defendiam ser uma decisão acertada, enquanto outros advogavam que os investimentos deveriam ser feitos em outras áreas da saúde.

Por curiosidade, acessei os dados da WHO (World Health Organization) e do banco mundial (World Bank) sobre quantidade de médicos por país e indicadores de saúde. Minha expectativa era encontrar pelo menos uma tímida relação entre indicadores. Mais do que isso, entender qual a localização do Brasil em relação a outros países. Fui surpreendido por uma forte correlação, que exploraremos a seguir.

Adotamos países como unidade observacional com medidas \(x\), o número de médicos 1,000 habitantes, e \(y\), a expectativa de vida saudável ao nascer.
Usando dados obtidos dos portais da WHO e do World Bank, plotamos os pontos no plano cartesiano.

    # http://apps.who.int/gho/data/view.main.HALEXv
    # https://data.worldbank.org/indicator/SH.MED.PHYS.ZS
    >library(magrittr)
    >library(ggplot2)
    >library(dplyr)    

    >worldbank_df <- read.csv("data/API_SH.MED.PHYS.ZS_DS2_en_csv_v2_10227587.csv",
                            header = T,skip = 3)
    >colnames(worldbank_df)[1] <- "Country"    

    >worldbank_df$n_docs <- sapply(split(worldbank_df[,53:62], #lists of values
                                        seq(nrow(worldbank_df))),
        function(x) tail(x[!is.na(x)],1)) %>% #ultimos valores não nulos
      as.numeric    

    >who_df <- read.csv("data/who_lifeexpect.csv",skip=2)
    >who_df$hale <- who_df$X2016
    >uni_df <- left_join(worldbank_df[,c("Country","n_docs")],
                        who_df[,c("Country","hale")],by="Country")    

    >ggplot(uni_df,aes(x=n_docs,y=hale))+
      geom_point(alpha=0.5,size=3) +
      xlab("No. de medicos / 1,000 hab.")+
      ylab("Expectativa de vida saudavel ao nascer")

É evidente que o padrão não é aleatório. Visualmente, notamos que o valor da expectativa de vida aumenta com maior Nº de médicos. Ainda, notamos um aumento inicialmente rápido até atingir um platô. O padrão é semelhante ao de uma curva logarítmica.

\(y = log(x)\) ou \(HALE = log(N_{médicos})\)

Se essa hipótese for verdade, transformar o número de médicos usando função logaritmica tornará a relação linear com a variável transformada:
Se \(y = log (x)\), fazemos a substituição \(x’ = log(x)\) para obtermos \(y = x'\).

Então a expectativa de vida se torna linearmente correlacionada ao logaritmo do número de médicos.

    >uni_df$log_docs <- log(uni_df$n_docs)
    >ggplot(uni_df,aes(x=log_docs,y=hale))+
      geom_point(alpha=0.5,size=3) +
      xlab("ln de No. de medicos / 1,000 hab.")+
      ylab("Expectativa de vida saudavel ao nascer")

De fato, verificamos uma notável tendência linear para os pontos.

Usando a implementação nativa em R para o coeficiente de Pearson:

    >cor.test(uni_df$log_docs,uni_df$hale)
    Pearson's product-moment correlation
    data:  uni_df$log_docs and uni_df$hale
    t = 18.572, df = 143, p-value < 2.2e-16
    alternative hypothesis: true correlation is not equal to 0
    95 percent confidence interval: 0.7854248 0.8828027
    sample estimates:
        cor
    0.8407869

A correlação linear obtida para nossa amostra de países é surpreendentemente grande, como sugeria a visualização \((\rho \sim 0.841)\).

O valor p é baixo \((p<0.001)\) considerando a hipótese nula \(H_{0}\) de \(\rho=0\). Concluímos então que há uma relação linear significativa de forte magnitude entre o logaritmo do número de médicos e a expectativa de vida dos países em nossa amostra.

É realmente curioso que exista uma relação matemática tão evidente entre construtos tenuamente conectados. O tempo médio que um organismo leva entre nascimento e morte e o número de profissionais atuantes. É virtualmente impossível explicitar cada relação causal por trás dessa relação, que se manifesta de forma robusta através da soma de muitos fatores relacionados.

Avaliando performance

Existem diferentes parâmetros para avaliar a performance de um modelo. Boa parte da pesquisa em aprendizagem estatística hoje consiste em implementar heurísticas que levem a melhores indicadores de performance.

Para regressão linear, o \(R^{2}\) é um coeficiente bastante usado. Expressa a proporção entre (1) variância explicada pelo modelo e (2) variação total. Chamamos de resíduo(ou erro) a diferença entre valores preditos e valores reais.

(1) Para capturar a magnitude dos erros do modelo, somamos o quadrado de todos os resíduos (sum of squared residuals, SSR) em relaçã aos valores preditos. Sejam \(y_{i}\) as observações e \(\hat{y_{i}}\) as predições:
\[SSR=\sum_{i=1}^{n} e_{i}^{2}=\sum_{i=1}^{n}(y_{i}-\hat{y_{i}})^{2}\]

(2) A variabilidade total é quantificada pela soma do quadrado dos desvios em relação à média (total sum of squares, TSS), um termo que vimos no cálculo da variância (segundo momento): \[TSS=\sum_{i=1}^{n}(y_{i}-\mu_{y})^{2}\]

Então a fração \(\frac{SSR}{TSS}\) é a proporção desejada. Definimos \(R^{2}\) como: \[R^{2}=1 - \frac{SSR}{TSS}\]

Uma visualização intuitiva de SSR e TSS:

    >source("aux/multiplot.R")
    >doc_lmfit <- lm(hale ~ log_docs, data=uni_df)
    >uni_df$preds[complete.cases(uni_df)] <- predict(doc_lmfit) 
    >uni_df$hale_mean <- mean(uni_df$hale,na.rm = T)
    >ssr_res <- ggplot(uni_df,aes(x=log_docs,y=hale))+
        geom_point(alpha=0.5,size=3) +
        geom_segment(aes(xend = log_docs, yend = preds)) +
        geom_smooth(method="lm")+
        xlab("")+
        ylab("Expectativa de vida saudavel ao nascer")+
        ggplot2::ggtitle("SSR")
    
    >tss_res <- ggplot(uni_df,aes(x=log_docs,y=hale))+
        geom_point(alpha=0.5,size=3) +
        geom_segment(aes(xend = log_docs, yend = hale_mean)) +
        geom_abline(slope = 0,intercept = 63.28165)+
        xlab("ln de No. de medicos / 1,000 hab.")+
        ylab("Expectativa de vida saudavel ao nascer")+
        ggplot2::ggtitle("TSS")
    
    >multiplot(ssr_res,tss_res)

O quadrado da distância entre um ponto e a reta corresponde a um resíduo. Obtemos SSR e TSS somando todos os resíduos nas figuras superior e inferior, respectivamente.

Valores de \(R^{2}\) próximos a 1 indicam soma de resíduos (SSR) similar a 0. Valores de \(R^{2}\) próximos a 0 indicam \(\frac{SSR}{TSS} \sim 1\) e as predições obtidas pelo modelo são tão boas quanto chutar a média para todos os casos.

    >lm(hale ~ log_docs, data=uni_df) %>% summary
    Call:
    lm(formula = hale ~ log_docs, data = uni_df)    

    Residuals:
         Min       1Q   Median       3Q      Max 
    -12.0964  -2.3988   0.3233   2.8229   8.6708     

    Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
    (Intercept)  64.4613     0.3162  203.84   <2e-16 ***
    log_docs      3.7303     0.2009   18.57   <2e-16 ***
    ---
    Signif. codes:  
    0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1    

    Residual standard error: 3.779 on 143 degrees of freedom
      (119 observations deleted due to missingness)
    Multiple R-squared:  0.7069,    Adjusted R-squared:  0.7049 
    F-statistic: 344.9 on 1 and 143 DF,  p-value: < 2.2e-16

Para obter os valores preditos, usamos o método predict:

    >head(predict(doc_lmfit))
        
           2        3        4        7        8        9 
    59.90747 57.23226 65.39962 66.11533 69.54483 68.30608 

É possível também obter predições para novos valores especificando o argumento newdata. Para um país com 1.5 médicos/1,000 habitantes:

    >predict(doc_lmfit,newdata = data.frame(log_docs=log(1.5)))
           1 
    65.97381 

Assumptions

Existem alguns procedimentos auxiliares para checar possíveis falhas e pontos no modelo que precisam de atenção. Por exemplo, os resíduos podem ser assimétricos. Isso indica que o desempenho muda em diferentes intervalos (heteroscedacidade). Uma lista completa de premissas, junto aos códigos em R para testá-las, está disponível no material auxiliar (lm-asssumptions.R)

Exercícios

1. Usando o dataset iris:

• Calcule medidas de tendência central e dispersão para cada variável.
  • Execute um teste t para uma das medidas entre duas espécies.
    
  • Obtenha o tamanho de efeito (D de Cohen) para a diferença.
    
• Faça um scatterplot entre duas medidas
  • Verifique se há correlação linear significativa entre as variáveis.
    
  • Se existir, ajuste um modelo de regressão linear.
    
  • Adicione cores de acordo com a espécie.
    
  • Ajuste um modelo de regressão para cada espécie.
    
  • Observe os valores de \(R^{2}\) para cada modelo. Qual a sua impressão sobre as mudanças de performance?

Ranks e o \(\rho\) Spearman

Relações lineares mantêm proporções constantes e aprendemos como quantificá-las. Por outro lado, duas variáveis podem ter relações de outros tipos, não lineares. Em especial, se as medidas apresentam valores muito extremos (outliers) um cálculo como o anterior sofre bastante com vieses.
Uma simples solução para esse problema é ranquear os valores. Assim, os itens do conjunto são tratados pela sua posição em relação a outros itens, de forma independente dos valores associados. Exemplo:
\[S = (1,3,89,89,39,209) \rightarrow S_{ranked} = (1,2,4,4,3,5)\]

O \(rho\) de Spearman é que o coeficiente produto-momento de Pearson aplicado aos ranks. Assim, medimos o grau em que duas variáveis aumentam (ou diminuem) em magnitude observando apenas a ordem das observações. Isto é: maior que, igual ou menor que. Especificamente, investigamos se há uma relação de monotonicidade entre elas.

Para a relação (sigmoide), entre x e y abaixo:

    >sig_data <- data.frame(y_vals = -(1 / (1 + exp(seq(-10,10,by =0.3) ) ) ),
                       x_vals = 1:67)
    >ggplot(sig_data,aes(x=x_vals,y=y_vals))+
    geom_point()

O coeficiente de Pearson é \(\rho \sim 0.936\)⁹ :

    >cor.test(sig_data$y_vals,
    +          sig_data$x_vals)    

        Pearson's product-moment correlation    

    data:  sig_data$y_vals and sig_data$x_vals
    t = 21.462, df = 65, p-value < 2.2e-16
    alternative hypothesis: true correlation is not equal to 0
    95 percent confidence interval:
     0.8978064 0.9603803
    sample estimates:
          cor 
    0.9361287  
    
    >ggplot(sig_data,aes(x=x_vals,y=y_vals))+
      geom_point()+
      geom_smooth(method="lm")

Como a relação é perfeitamente monotônica, os pares ordenados \((x_{i},y_{i})\) sempre possuem o mesmo rank. O quinto valor mais alto em x é também o quinto valor mais alto em y. Portanto, o coeficiente de Spearman é 1:

    >cor.test(sig_data$y_vals,
    +          sig_data$x_vals,method = "spearman")    

        Spearman's rank correlation rho    

    data:  sig_data$y_vals and sig_data$x_vals
    S = 0, p-value < 2.2e-16
    alternative hypothesis: true rho is not equal to 0
    sample estimates:
    rho 
      1 

O coeficiente \(\rho\) de Spearman é preferível quando as medidas diferem bastante de uma distribuição normal. Também oferece robustez quando é necessário lidar com outliers.